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Abstract 

A free-energy functional that contains both the symmetry conserved and symmetry broken parts 
of the direct pair correlation function has been used to investigate the freezing of a system of 
hard spheres into crystalline and amorphous structures. The freezing parameters for fluid-crystal 
transition have been found to be in very good agreement with the results found from simulations. 
We considered amorphous structures found from the molecular dynamics simulations at packing 
fractions n lower than the glass close packing fraction nj and investigated their stability compared 
to that of a homogeneous fluid. The existence of free-energy minimum corresponding to a density 
distribution of overlapping Gaussians centered around an amorphous lattice depicts the deeply 
supercooled state with a heterogeneous density profile. 

PACS numbers: 64.70.D-, 05.7Q.fh, 64.70.pm 
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I. INTRODUCTION 



When a liquid freezes into a crystalline solid its continuous symmetry of translation and 
rotation is broken into one of the Bravais lattices. In three dimensions the freezing of a liquid 
into a crystalline solid is known to be first-order symmetry breaking transition marked by 
large discontinuities in entropy, density and order parameters. The molecules in a crystal are 
localized on a lattice described by a discrete set of vectors {Rj} such that any functions of 
position, such as one particle density p(r) satisfies p(r) = p(r + Rj) for all Rj [1]. This set of 
vectors necessarily forms a Bravais lattice. But when the liquid is supercooled bypassing its 
crystallization, it continues to remain in amorphous state. With increase in density a solid 
like phase emerges with molecules getting localized around their mean position on a random 
structure. The underlying lattice on which such localized motion takes place is related to 
the time scale of relaxation in the supercooled liquid. While the supercooled liquid starts 
to attain solidlike properties, structurally it does not have any long range order like the one 
present in the crystal. This phenomena is termed the "glass transition" [2,3]. Although 
the glassy materials are well characterized experimentally, the existence of phase transition 
into the glassy state remains controversial [4 — 6]. Our aim in this paper is not to enter in 
this discussion but to examine the stability (metastability) of amorphous structures from a 
thermodynamic point of view, using the standard method of density-functional theory which 
is also used to investigate the crystallization of liquids. We refer to glassy or amorphous solid 
a structure in which particles are localized around their mean positions forming a random 
lattice. Localization of particles amounts to breaking of continuous translational symmetry 
of normal liquid and takes place in forming both crystals and glasses. 

The structural properties of a classical system can adequately be described by one and 
two-particle density distributions. The one particle density distribution p(r) defined as 



where Rk is the position vector of k-th particle and bracket (....) represents the ensemble 
average, is constant independent of position for an isotropic liquid, but contains most of the 
structural informations of crystals and glasses. The particle density distribution p( 2 )(r!,r 2 ) 
which gives the probability of finding simultaneously a molecule in a volume element dr\ 
centered at ri and a second molecule in volume element dr 2 centered at r 2 is defined as 



p(r) = (£<Hr-R fe )> 




k 



2 



(ri , r2 ) = 5( Tl - Rj)5(r 2 - R k )> (1.2) 

3 

The pair correlation function (?(ri,r 2 ) is related to p( 2 )(ri,r 2 ) by the relation 

p( 2 )( ri ,r 2 ) 

#ri,r 2 = — — — — 1.3 
P(ri)p(r 2 ) 

The direct pair correlation function c(ri, r 2 ) which appears in the expression of free-energy 
functional (see Sec. II) is related to the total pair correlation function h(ri, r 2 ) = g(r 1 , r 2 ) — 1 
through the Ornstein-Zernike (OZ) equation, 

/i(ri, r 2 ) = c(ri, r 2 ) + J c(r u r 3 )p(r 3 )/i(r 2 , r 3 )dr 3 (1.4) 

Since in an isotropic liquid p(ri) = p(y 2 ) — Pi — N/V where N is the average number of 
molecules in volume V, 

p (2) (n,r 2 ) 

#( r ) = — m — ( L5 ) 

Pi 

where r = |r 2 — r i I • 111 a liquid of spherically symmetric particles g{r 1 ,v 2 ), h(r 1 ,r 2 ), c(r 1? r 2 ) 
depend only on interparticle separation r = |r 2 — ri|. This simplification is due to homogene- 
ity which implies continuous translational symmetry and isotropy which implies continuous 
rotational symmetry. Such simplification does not, in general, occur in frozen phases. We 
refer crystal as well as glass as frozen phases. While a crystal is both inhomogeneous and 
anisotropic, a glass can be regarded as isotropic but inhomogeneous. The heterogeneity in 
glassy system over length and time scales has been studied in several recent works [7] related 
to computer simulations. 

The total and direct pair correlation functions of a system can be given as a simultaneous 
solution of the OZ equation and a closure relation that relates the correlation functions to 
the pair potential. Well known closure relations are the Percus-Yevick (PY) relation, the 
hypernatted chain (HNC) relation and the mean spherical approximation (MSA). It may, 
however, be noted that while the OZ equation (1.4) is general and connects the total and 
direct pair correlation functions of liquids as well as of frozen phases, the closure relations 
have been derived assuming translational invariance [8] and therefore valid only to normal 
fluids. The integral equation theory has been used quite successfully to describe the structure 
of isotropic liquids. But its application to frozen phases has so far been very limited [9, 10]. In 
Sec. Ill we describe a method to calculate the DPCF of frozen phases formed by breaking of 



continuous translational symmetry of liquids and use it in a free-energy functional described 
in Sec. II to study freezing transitions in Sec. IV and V. 

Since its advent in 1979 by Ramakrishnan and Yussouff (RY) [11], the density functional 
theory (DFT) has been applied to freezing transition of variety of pure liquids and mixtures 
[11, 12]. A DFT requires an expression of the Halmholtz free-energy (or grand thermody- 
namic potential) in terms of one- and two-particle distribution functions and a relation that 
relates p(r) to pair correlation functions. Such a relation is formed by minimizing the free- 
energy with respect to p(r) with appropriate constraints [12, 13]. The DPCF that appear 
in these equations are of frozen phase and are functional of p(r). When this functional 
dependence is ignored and the DPCF is replaced by that of the co-existing isotropic liquid 
[11] or by an "effective fluid" [14] the free-energy functional becomes approximate and fails 
to provide an accurate description of freezing transitions and stability of frozen phases. An 
improved free-energy functional which contains both the symmetry conserved and symme- 
try broken parts of the DPCF has recently been developed [9, 10] and applied to study the 
isotropic-nematic transition [9] and the crystallization of power-law fluids [10]. 

In this paper we investigate the freezing of fluids of hard sphere into crystalline and 
amorphous phases and compare our results with the results of previous investigations. Hard 
spheres are ubiquitous in condensed matter; they have been used as models of liquids, 
crystals, colloidal systems, granular systems and powders. Packing of hard spheres are of 
even wider interest as they are related to important problems in information theory, such as 
diagonalization of signals, error correcting codes, and optimization problems [15]. Recently, 
amorphous packing of hard spheres have attracted much attention [6, 16 — 18] because for 
polydisperse colloids and granular materials the crystalline state is not obtained for kinetic 
reasons . It is therefore necessary to have a statistical-mechanical theory based on first 
principle which can correctly describe the freezing of hard spheres. 

The paper is organized as follows: In Sec. 77 we describe the free-energy functional for a 
symmetry broken phase which contains both the symmetry conserving and symmetry broken 
parts of direct pair correlation functions. In Sec. Ill we describe a method to calculate 
these correlation functions. The theory is applied in Sec. IV to investigate the freezing of 
fluids into crystalline solids and in Sec. V the metastability of amorphous structures. 
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II. FREE-ENERGY FUNCTIONAL 



An important step in construction of the density functional model of a frozen phase is 
proper parametrization of extremely inhomogeneous density function p(r) which value near 
a lattice site may be orders of magnitude higher than in the interstitial regions. One very 
successful prescription of p(r) is as a collection of overlapping Gaussian profiles [19] centered 
over a set of lattice sites R m . 



P(r) = E(-) 3/2 e^ -a(r - R m ) 2 (2.1) 



m 



3/2 ( 

IT' 



The width of a Gaussian profile is inversely proportional to the square root of a which 
will be referred to as localization parameter. In this representation, the limit a — > is the 
homogeneous liquid state and higher values of a represent increasingly localized structures. 

Taking Fourier transform of (2.1) one gets 

P(r)=Po + ^E/V^ r (2.2) 

where 

Pg = e - 9 2 /4a^ e -,qR„ (2.3) 

n 

is amplitude of density wave of wavelength 2ir/\q\. The nature and magnitude of inhomo- 
geneity of a frozen phase is measured by p q which will be referred to as order parameter; 
p q = for q 7^ corresponds to isotropic fluid and p q ^ to a frozen phase. For a crystal 
in which R m forms a periodic lattice, e* q Rm = 5^g where G are reciprocal lattice vectors 
(R. L. V.), Eq.(2.2) reduces to 

p s (r)=p + p Y / e' G2/4a e^ (2.4) 

G 

This is a well known expression of p(r) of a crystal. As in an amorphous or glassy structure 
the lattice sites are randomly distributed and are not known, the above simplification is not 
possible. We can, however, calculate and have an idea of inhomogeneity and its difference 
from that of a crystalline solid using the distribution of R n determined from numerical 
simulations [20,21]. Thus 



P 9 (r) = Po + ^ E e ^ 2/4 ° \ l + Po ( I dRh(R)e-^ 



e jqr (2.5) 
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where g(R) = 1 + h(R) is the site-site correlation function which provides the structural 
description of amorphous structure. In Fig.l we plot p g (r) for r](= ^npoa 3 ; a being the 
diameter particle)= 0.576 and a = 150 (highly localized condition) and a — 15 (weakly 
localized condition). In calculating p g (r) we used g(R) data found for amorphous structures 
of granular particles [21] shown in Fig.2. In Fig. 2 we also show the value of g(R) of a liquid 
found from solving the integral equations discussed in Sec. Ill, for comparison. From Fig.l 
we see that while inhomogeneity increases with value of a, unlike crystalline solid, it remains 
confined to about 4 particle diameter, even for a = 150 and decays rapidly on increasing 
the distance. 

The reduced free-energy A[p] of an inhomogeneous system is a functional of p(r) and is 
written as 



A[p] = A id [p] + A ex [p] (2.6) 
The ideal gas part is exactly known and is given as 

A id [p\- / drp(r) Mp(r)A) - 1] (2.7) 

where A is cube of the thermal wavelength associated with a molecule. The excess part 
arising due to intermolecular interactions is related to the DPCF as [9, 10] 

_|^d_._ c (o» (ri , r2;po) _ c ,, (ri , r2;[p]) (2 . 8) 

where superscripts (0) and (6) represent, respectively, the symmetry conserving and sym- 
metry breaking parts of the DPCF. In other words is found by treating the system to 
be isotropic and homogeneous with density po whereas are the contribution which arise 
due to heterogeneity in density in a frozen phase. 

A e;E [p] is found by functional integration of Eq.(2.8). In this integration the system is 
taken from some initial density to the final density p(r) along a path in the density space; 
the result is independent of the path of integration. As the symmetry conserving part 
is function of density the integration in the density space is done taking an isotropic fluid of 
density p (or pi, the density of co-existing fluid in case of crystal) as reference. This leads to 

45 M = A ex (p ) rfr 2 Ap(r 1 )Ap(r 2 )c(°)(r) (2.9) 
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where Ap(ri) = p(ri) — po an d A ex (p ) is the excess reduced free-energy of an isotropic 
system of density po- 

The integration over has to be done in the density space spanned by the number 
density po and order parameters p q . We characterize this part of the density space by two 
parameters A and £ which vary from to 1. The parameter A raises the density from to 
Po as it varies from to 1 whereas the parameter £ raises the order parameters from to p q 
for each value of q. This gives [9, 10]. 

42 M = ~\j J dr 2 Ap( ri )Ap(r 2 )cW (n,r 2 ; [p]) (2.10) 



where 

7&)t 



•!,r 2 ;[p]) = 4 fdSS /V f 1 dXX C dX'c^, r 2 ; AA'p ; tf' p g ) (2.11) 

JO JO JO JO 

While integrating over A the order parameter p q are kept fixed and while integrating 
over £ the density is kept fixed. The result does not depend on the order of integration. 
The free-energy functional of a frozen phase is the sum of A id , AA^ [p] and AA^ [p] given, 
respectively, by Eq. (2.7), (2.9) and (2.10). 

The minimization of AA = A[p] — A(p$) where A(po) is the free-energy of an homogeneous 
and isotropic system of density po, leads to 

pW 

Po 

Here A is the Lagrange multiplier and is determined from the condition 

1 f P{r) -dr=l (2.13) 



ln^£) =A + / rfr 2 Ap(r 2 )c(°)(|r 2 -r 1 |,p )+ / rfr 2 Ap(r 2 )g( b )( ri , r 2 ) (2.12) 
On J J 



V Jv po 
and 

c^( ri ,r 2 ) = 2 TdX f 1 d^ir^Xp^M (2.14) 
JO jo 

In principle, the only information we need to know is the value of c^(r) and c^(ri, r 2 ; \p\) 
to calculate self consistently the value of p(r) which minimizes the free-energy. In practice, 
one, however, finds it convenient to do minimization for assumed form of p(r) [12] 

III. CALCULATION OF DIRECT PAIR CORRELATION FUNCTION 

To calculate the values of c^(r) we use the integral equation theory consisting of OZ 
equation, 

/i(°)(r) = c (°)(r)+p / drV 0) (r>^(|r'-r|) (3.1) 



and a closure relation proposed by Rogers and Young [22]. This closure relation is written 



as 



1 + h {0 \r) = exp{-u{r)/k B T) 

where 



l exp{l{r)f{r)) 



(3.2) 



fir) 

7 (r) = h(r) - c(r) (3.3) 



and 

/(r) = 1 — exp(-ipr) (3.4) 

Here ^ is an adjustable parameter used to achieve thermodynamic consistency and its 
value for a system of hard spheres is found to be equal to 0.16 [22]. In Eq.(3.2) u(r) is 
pair potential, k B , the Boltzmann constant and T, temperature. This closure relation was 
found by mixing the PY and HNC closure relations in such a way that at r = it reduces 
to the PY approximation and for values of r where /(r) approaches 1, it reduces to the 
HNC approximation. Eqs(3.1) — (3.4) together constitute a thermodynamically consistent 
theory and has been found to give values of pair correlation functions which are in very good 
agreement with Monte Carlo results. 

The differentiation of Eqs.(3.1) and (3.2) with respect to density yields the following two 
relations, 

/ rfr'c(°)(r> (|r'-r|)+po / dr'^h^\\r' -r\)+p f dr'c<°>(r) ~ ^ 

(3.5) 



dp dp 
and 



+ 



dh^ir) ( u{r) 



exp 



exp[~f(r)f(r)] 



dpc 



(3.6) 



dp V k ^T / 

The closed set of coupled equations (3.1), (3.2), (3.5) and (3.6) have been solved using 
Gillen's algorithm [23] for four unknowns h^°\ c^°\ and We compare the values 

of c(°)(r) at packing fractions rj = 0.50 and 0.55 in Fig. (3) and values of 8c ^ r - > for the same 
two values of rj in Fig. (4) to see the density dependence of these functions. While rj = 0.50 is 
close to the value of packing fraction at which system freezes into a crystalline solid, rj = 0.55 
is close to the value of packing fraction t)m at which the crystal melts. 

Since all closure relations (including the one given by Eq.(3.2)) which are used in the 
integral equation theory for pair correlation functions are derived assuming translational 
invariance [8], their use in calculating the values of pair correlation function of frozen phases 
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may not be appropriate. In view of this we use a series expansion in which the contribution 
to the the DPCF arising due to inhomogeneity of the system is expressed in terms of higher 
body direct correlation functions of the uniform (isotropic and homogeneous) system. Thus 

c (b) (r!,r 2 ; [p]) = J dr 3 4 0) (n, r 2 , r 3 ; p )(p(r 3 )-p ) + J dr 3 dr 4 c{ 0) (r u r 2 , r 3 , r 4 )(p(r 3 )-p )(p(r 4 )-p )+- 

(3.7) 

where p(r n ) — p = J2 q ^o P 9 e iqrn . In Eq.(3.7) are n-body direct correlation functions of 
the uniform system which can be found using the relations 



Spo 



5 2 c (o )( 



r 



= J A 3 cf (n,r 2 ,r 3 ) (3.8) 

A2 - /dr 3<i r 1 cr(r 1 ,r 2 ,r„r 4 ) (3.9) 
dp J 

etc. The values of derivatives of c^°^(r) appearing on 1. h. s. of above equations have 
been found using the integral equation theory described above. 

We note that Eq.(3.7) satisfies the condition that is zero in the fluid phase and 
depends on the magnitude (order parameter) and phase factors of density waves. These 
density waves measure the nature and magnitude of inhomogeneity of frozen phases. While 
each wave contributes independently to the first term of Eq.(3.7) interactions between the 
two waves contribute to the second term and so on. The contributions made by successive 
terms of Eq.(3.7) depends on the range of pair potential u(r) [10]. As the range of potential 
increases the contribution made by higher order terms increases. For a system of hard 
spheres we find that at the freezing transition the contributions made by first term to free- 
energy is already small and therefore higher terms are expected to be negligible; it is only 
for u(r) oc r~ n ,n < 12 the contribution made by second order term becomes important 
[10]. In view of fast convergence of the series, Eq.(3.7) seems to be a useful expression for 
calculating c^(i"i,r 2 ) 

The first term of Eq.(3.7) involves three-body direct correlation function which can be 
factorized as a product of two-body functions [24] . Thus 

4 0) (ri,r 2 ,r 3 ;A>) = ^12)^13)^23) (3.10) 
The function t(r) is determined using relation of Eq. (3.8) 

?^l = t(r)Jdr't(r)t(\r-r\) (3.11) 



We adopt the numerical procedure developed in [24] to calculate t(r) from known values 
of ^j^p- from (3.11). The values of t(r) are plotted in Fig. (5) for rj = 0.50 and 0.55 to show 
their density dependence. 

Taking only the first term in Eq.(3.7) we write, 



(ri,r 2 ) = ^EE^/ drstdrs-rxDe^-^drs-ral) (3.12) 



where [i q = e q2 ^ 4a Using the relation 

r 3 = ^(ri + r 2 ) - ^(r 2 - r x ) + (r 3 - r x ) (3.13) 

we find that c^(ri,r 2 ) can be written in a Fourier series in the center of mass variable 
r c = ||ri + r 2 | with coefficients that are functions of the difference variable r = r 2 — ri, i.e, 

c (b) (ri,r 2 ) = ^E c(9) We iqrc (3.14) 

where 

c («)(r) = EA i ? e ^ q ' Rne_ ^ qr / rfr't(r')e iq - r 't(|r' -r|) (3.15) 

Since the DPCF is real and symmetric with respect to the interchange of i*i and r 2 , 
c^(r) = c(~ q \r) and c^(r) = c^(— r). For given value of a and R n one can calculate 
c (?)( r ) from known value of t(r). We discuss our results for crystalline and amorphous solids 
in following sections. 

IV. CRYSTALLINE SOLID 
A. Calculation of c (b) (r!,r 2 ) 



For a crystal in which vectors R n form a Bravais lattice Eq. (3.14) and (3.15) can be 
written as [10] 

c ( fe )( ri ,r 2 ) = Ec (G) (r)e lG ' rc (4.1) 

G 

and 

c ( G )(r) = p ii G e-t iG r J dr't(r')e iG - r 't(\r - r|) (4.2) 
where /iq = e~ G2 ^ 4a Eg. (4.2) has been solved using Rayleigh expansion to give 

c (G) (r)=Y,€\r)Y lm (r) (4.3) 

Ira 
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where 



r) 



2tt 2 



EEO 

h Z 2 



-11 



(2/ 1 + l)(2/ 2 
2/ + 1 



[C g (h, l 2 , 1; 0, 0, 0)] 2 J/2 (-Gr)t(r)B h (r, G)^(G) 



(4.4) 



Here C 9 is the Clebsch-Gordan coefficient, j t (x) the spherical Bessel function and 



B h (r,G) = (An) 2 J dkk 2 t(k)j h (kr) J dr ' r 2 t(r') 3h (kr) 3h (Gr') (4.5) 

The crystal symmetry dictates that / and h+h are even and for a cubic crystal, m = 0, ±4. 
The c^(r) depends on the order parameter and on the magnitude of R. L. V. 
The Fourier transform of (r) defined as 



c<°>(k) = J c^\v)e-^dv = J2€\WUk) 

lm 



(4.6) 



where 



c\Z\k) = J drr 2 3l (kr)c^\r) (4.7) 

is calculated from the knowledge of c lm (r). In Figs(6) we plot c lm (k) for the first two sets 
of G at rj = 0.55 and a = 170 for face centered cubic (f.c.c.) structure. 
B. Liquid-Solid transition 

The grand thermodynamic potential defined as —W — A — /3/j, J drp(r), where fx is the 
chemical potential, is used to locate transition as it ensures that the pressure and chemical 
potentials of two phases remain equal at the transition. The transition point is determined 
by the condition AW = Wi — W = where W\ is the grand thermodynamic potential of the 
fluid. Using expressions given in Sec. II we find 



AW AW id AWo AW b 



N 



N 



N 



N 



(4.8) 



where 



AW, 



id 



N 



= l-lnpi + (l + Ap* 



3, 5 

2 in y-2 



(4.9) 



N 

AWb 

N 



E'E'(! + Ap*) Wc-c^Gi + 



Gi G 



(4.10) 
(4.11) 
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Here AWid, AWq and AWb are, respectively, the ideal, symmetry conserving and sym- 
metry broken contribution to AW, the prime on summations in Eq.(4.11) indicates the 
condition, G ^ 0, Gi ^ and G + G x ^ 0, and 

c(°)(G) = Pl J c {Q \r)e tG - r dv (4.12) 
S(d + ±G) = po E / PoK (Gl ^ G > r >W^r (4.13) 

lm 

where p = pi(l + Ap*) 

We used the above expression to locate the fluid-f.c.c. solid and fluid-b.c.c. solid transi- 
tions by varying the values of pi, Ap* and a. While we find fluid-f.c.c. solid transition to 
take place at rji = 0.490, Arj* = 0.106 and a = 170, no transition is found for b.c.c. solid. 
In table 1 we compare our results of freezing parameters with those found by Monte Carlo 
simulation [25,26] and from other density functional schemes [27-30]. The agreement found 
between our results and those of simulations are very good, better than any other density 
functional schemes. 

At the transition point the contribution of different terms of Eq.(4.8) is as follows; AWa — 
4.44, AWo — —4.10, AWb _ _o.34. The contribution made by symmetry breaking term of 
free-energy is about 8.3% to that of the symmetry conserving term. This is in accordance 
with the result found earlier [10] for inverse power potential u(r) = e(a/r) n , where e, a 
and n are potential parameters, that as n increases(n = oo corresponds to hard sphere 
potential) the contribution made by the symmetry breaking term to free-energy decreases. 
This explains why RY theory [11] while gave good results for hard spheres system, failed for 
potentials which have soft repulsion and/or attractive tails. 



V. AMORPHOUS STRUCTURE 



In this section we investigate the heterogeneous density profile of an amorphous structure 
and examine the question of having metastable states in between the normal fluid state and 
the regular crystalline state at packing fraction r\ which lies between packing fraction at 
melting point of a crystal r\ M = 0.545 and packing fraction corresponding to "glass close 
packing", r]j ~= 0.65. The usual way to construct amorphous structure in experiment or 
numerical simulation is to compress the system according to some protocol which avoids 
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crystallization [17,20,21]. One of the criteria used to signal the onset of glassy phase in 
supercooled liquids is emergence of split second peak. There may be infinitely large number 
of such metastable structures which when compressed jam along a continuous finite range 
densities down to the glass close packing rjj [31, 32]. 

The density functional approach provides the means to test if such a structure is stable 
compared to that of a fluid at a given temperature and density. In earlier calculations the 
random closed packed structure generated through Bennett's algorithm [33] was used. The 
g(R) giving the distribution of particles at a given value of r\ was found using an ad hoc 
scaling relation [34]. 

g(R)=g B [R(^] (5.1) 

where r/j was used as a scaling parameter such that at i] = r/j the random closed packed 
structure <?_b(-R) was obtained. While Singh et. al. [35] found that the state corresponding to 
this structure becomes more stable than fluid for rj > 0.59 for very large value of a (~ 280) 
that corresponds to highly inhomogeneous density distribution, Kaur and Das [36] found that 
the same structure also becomes more stable than fluid for rj > 0.576 for considerably smaller 
values of a (~ 18). On the other hand, Dasgupta [37] has numerically located the "glassy" 
minimum of a free-energy functional and the structure which gave this minimum. Here 
we use the value of g(R) found for granular particles from molecular dynamics simulations 
[21] at rj = 0.576 and 0.596 and examine the stability of these structures. The reason for 
choosing these data is that they are available at rj < rjj and have the essential features of 
an amorphous structure. 

A. Calculation of c( b )(n,r 2 ) 

From the known values of g(R) the order parameter defined by Eq.(2.3) is calculated. 
Thus, for q ^ 

Pq = ^ 2 ^e~^= N S a (q) (5.2) 

n 

where n q = e - q2 ' ia and S q (a) = 1 + 24r? / dRR 2 (g(R) - l)j (qR). In Fig. (7) we plot p q as 
a function of q for a = 15 and 50 and i] = 0.596. The value of order parameter /i G of a 
f.c.c crystal for a = 50 and i] = 0.596 are also shown for comparison. From the figure one 
may note that the value of p q of an amorphous structure has a very different magnitude and 
dependence on wave vector q than that of a crystal. 
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Using the fact that an amorphous structure can be considered on the average to be 
isotropic, Eq.(3.15) is simplified to give 

c {q \r) = ^MtiEW + mMrtlqr) (5.3) 

In Fig. (8) we show the value of c^ q \r) as a function of r for q at which p q S a (q) is maximum 
for a = 15 and rj = 0.596. 

B. Determination of free-energy minimum 

We calculate the minimum of AA[p] = A[p] — Ai[pi] where Ai(pi) is the reduced free- 
energy of an isotropic fluid of density p and A[p] is the reduced free-energy of an amorphous 
structure of average density po- Using expression given in Sec. II we get 

AA[p] _ AAM AAM AA b [N] 
N N N N 1 ' ' 



^ =^{-)' 2 I rfrrV^ 2 ln[{(^) 3/2 e— 2 }+^(^) V2 [ dRRg(R){e'^ 2 -e^*) 2 }] 
N ix J ix r i J 

(5.5) 

for a < 20 

^ + -f] (5.6) 

for a > 20 

^ = ~i:\N\ 2 Sa(q)^(q) (5.7) 
= -^EE/W-^A(Si)c (9) (|qi + ^ql) (5-8) 



where 



and 



; qi + I q ) = J c {q \r)e l ^ + ^ r dr 



cW(r) =ljd\\j dX'dtf J dS'cto\r,\\'po,tt'p q ) 
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As shown in Fig(9), a minimum of free-energy is found at a ~ 7 for both r\ = 0.576 and at 
rj = 0.596. The width of a Gaussian density profile 1/y/a ~ 0.37. This is a case of weak 
localization. From the figure it is clear that the "glassy minimum" is separated from the 
liquid minimum by a barrier located at a ~ 0.6 which height grows with the density. The 
free-energy minimum corresponds to a density distribution of overlapping Gaussians centered 
around an amorphous lattice and depicts the deeply supercooled state with heterogeneous 
density profile. 

The schematic phase diagram that one expects in the presence of a glass transition, 
contains a pressure line that bifurcates from that of the liquid at some rj > tjm and which 
diverges at r]j [6]. The bifurcation point is connected with the onset of glass transition. The 
pressure of a system can be found from the knowledge of single and two particle density 
distributions. For example, the virial pressure of a system in 3-dimensions is given as 

(3P Iff 

= 1 - - / dr 1 / dr 2 p(ri)p(r 2 )0(ri, r 2 )(r • \ju{r)) (5.9) 

where r = r 2 — r 1 . For an amorphous structure of hard spheres this reduces to 

— = 1 + \ W(l) + ^(1) E \N\ 2 Sa(q)jo(q) (5.10) 

Note that for isotropic fluid the third term of above equation is zero; the bifurcation of 
pressure line from that of normal fluid starts as soon as particles start getting localized. 
Localization of particles also leads to crossover from nonactivated to activated dynamics 
and considerable increase in relaxation time. 



VI. SUMMARY AND PERSPECTIVES 



A free-energy functional that contains both the symmetry conserved part of the DPCF 
c(°)(r) and the symmetry broken part c^(r 1; r 2 ) has been used to investigate the freezing 
of a system of hard spheres into crystalline and amorphous structures. The values of (r) 
and its derivatives with respect to density p as a function of distance r have been found 
using integral equation theory comprising the OZ equation and a closure relation proposed 
by Roger and Young [22]. For c^(ri, r 2 ) we used an expansion in ascending powers of order 
parameters. This expansion involves higher-body direct correlation functions of the isotropic 
phase which in turn was found from the density derivatives of c^(r). For this we used the 

15 



ansatz [24] embodied in Eqs.(3.10) and (3.11). The contribution made by the symmetry 
broken term to the free-energy at the freezing point (liquid-crystal transition point) was 
found to be about 8% of the symmetry conserving part. Though this contribution is small 
but, as shown in Table 1, it improves the agreement between theoretical values of the freezing 
parameters and the values found from simulations. This result and the results reported 
earlier [10] for the power law fluids show that the contributions of the symmetry broken 
part of free-energy increases with the softness of the potential. This explains why the RY 
free-energy functional was found to give reasonably good description of freezing transition 
of hard spheres fluid but failed for potentials which have soft core and/or attractive tail. 
These results also indicate that the theory described here can be used to describe the freezing 
transitions of all kind of potentials. 

We used the free-energy functional to investigate the question of having metastable states 
in between the homogeneous liquid and the regular crystalline state. The value of site- 
site correlation function g(R) which provides the structural description of the amorphous 
structure have been taken from molecular dynamics simulation of granular system subjected 
to a sequence of vertical tapes [21]. The system has been found to behave like a glass- forming 
system. The reason for our choosing this data is that they are available at r) < rjj [21] and 
therefore can be directly used in the theory without using approximations such as scaling 
relation (5.1). Using the data of g(R) at rj = 0.576 and 0.596 from Ref [21] we examined 
the stability of amorphous structures with respect to homogeneous fluid. The minimum of 
free-energy found at a ~ 7 suggests that the structures are stable compared to that of the 
fluid and corresponds to a density distribution of overlapping Gaussians centered around an 
amorphous lattice. This kind of structure may be associated with deeply supercooled states 
with a heterogeneous density profile. The transition of the liquid into any of these states 
will be determined by considering the dynamics of fluctuations around these minima. The 
glassy minimum is separated from the homogeneous liquid minimum by an energy barrier 
which height increases with the density. 

Among future applications it will be instructive to investigate the contribution made 
by the second term of Eq.(3.7) to free-energy of different potentials, application of the 
theory to freezing transition in two dimensions, in particular to examine the melting through 
hexatic phase which other density functional schemes have failed to show and to examine 
the possibility of calculating total pair correlation functions of frozen phases using the OZ 
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equation. 

Acknowledgments: We are thankful to Dr. Massimo Pica Ciamarra and Dr. A. Donev 
for providing the data of g(R). One of us (S. L. S.) is thankful to the University Grants 
Commission for research fellowship. 



[1] P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics (Cambridge 

University Press, 1995) 
[2] M. D. Ediger, C. A. Angell and S. Nagel, J. Phys. Chem. 100, 13200 (1996) 
[3] P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001) 

[4] L. O. Hedges, R. L. Jack, J. P. Garrahan and D. Chandler, Science 323, 1309 (2009) 

[5] A. Cavagna, Phys. Rep. 476, 51 (2009) 

[6] G. Parisi and F. Zamponi, Rev. Mod. Phys. 82, 789(2010) 

[7] W. Kob., C. Donati, S. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett.79, 
2827(1997) 

[8] J. P. Hensen and I. R. Mc Donald, Theory of Simple Liquids, 3rd edition (Elsevier, 2006) 

[9] P. Mishra and Y. Singh, Phys. Rev. Lett. 97, 177801 (2006) 
[10] S. L. Singh and Y. Singh, Europhys. Lett. 88, 16005 (2009) 
[11] T. V. Ramakrishanan and M. Yussouff, Phys. Rev. B. 19, 2775 (1979) 
[12] Y. Singh, Phys. Rep. 207, 351 (1991) 
[13] H. Lowen, Phys. Rep. 237, 249(1994) 

[14] A. R. Denton and N. W. Ashcroft, Phys. Rev. A 39, 4701 (1989); A. Khen and N. W. Ashcroft, 

Phys. Rev. Lett. 78, 3346 (1997) 
[15] J. H. Conway and N. J. A. Sloane, Sphere packings, Lattices and Groups (Springer- Verlag, 

New York, 1993) 

[16] S. Torquato, T. M. Truskett and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000) 

[17] C. S. O'Hern, S. A. Langer, A. J. Liu and S. Nagel, Phys. Rev. Lett, 88, 075507 (2002); L. E. 

Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E. 73, 041304(2006) 
[18] R. D. Kamien, A. J. Liu, Phys. Rev. Lett. 99, 155501 (2007) 
[19] P. Tarazona, Mol. Phys. 52, 81 (1984) 

[20] A. Donev, F. H. Stilinger and S. Torquato, Phys. Rev. Lett. 96, 225502 (2006), A. Donev, R. 

17 



Connelly, F. H. Stillinger and S. Torquato, Phys. Rev. E. 75, 051304 (2007) 
[21] M. Pica. Ciamarra, M. Nicodemi and A. Coniglio, Phys. Rev. E 75, 021303(2007) 
[22] F. J. Rogers and D. A. Young, Phys. Rev. A. 30, 999 (1984) 
[23] M. J. Gillan, Mol. Phys. 38, 1781(1979) 

[24] J. L. Barrat, J. P. Hansen and G. Pastore, Mol. Phys., 63, 747(1988), Phys. Rev. Lett. 58, 
2075(1987) 

[25] W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968) 

[26] B. J. Alder, W. G. Hoover and D. A. Young, J. Chem. Phys. 49, 3688(1968) 

[27] D. C. Wang and A. P. Gast, J. Chem. Phys. 110, 2522(1999) 

[28] B. B. Laird, and D. M. Kroll, Phys. Rev. A. 42, 4810(1990) 

[29] J. L. Barrat, J. P. Hansen and G. Pastore and E. M. Waisman, J. Chem. Phys. 86, 6360(1987) 

[30] B. B. Laird, J. D. McCoy and A. D. J. Heymat, J. Chem. Phys. 87, 5449(1987) 

[31] L. Berthier and T. A. Witten, Phys. Rev. E. 80, 021502(2009), Y. Burmer, D. R. Reichmann, 

Phys. Rev. E. 69, 041202(2004) 
[32] P. Chaudhari, L. Berthier and S. Sastry, Phys. Rev. Lett. 104, 165701(2010) 
[33] C. Bennet, J. Appl. Phys. 43, 2727(1972) 

[34] M. Baus and Jean-Louis Colot, J. Phys. C 19, L135(1986) H. Lowen, J. Phys. Cond Matt 2, 
8477(1990) 

[35] Y. Singh, J. P. Stossel and P. G. Wolyness, Phys. Rev. Lett. 54, 1059(1985) 
[36] C. Kaur and S. P. Das, Phys. Rev. Lett. 86, 2062(2001) 

[37] C. Dasgupta, Europhys. Lett. 20, 131(1992); C. Dasgupta and O. T. Vails, Phys. Rev. E. 59, 
3123(1999) 



18 



TABLE I: Freezing parameters of a hard-sphere fluid derived from the various density functional 
schemes. Here L = ^ (l^r) ^ * s the Lindemann parameter, r] s = ^p s 0" 3 an d r\i = ^pia 3 . 
Average errors are given in the parentheses. MWDA stands for modified weighted density approx- 
imation and RY DFT stands for Ramakrishnan-Yussouff density functional theory. 







m 


At]* 


L 


Present result 


0.542(< 1%) 


0.490(< 1%) 


0.106(< 2%) 


0.09 


MWDA-static reference [27] 


0.503(8%) 


0.452(8%) 


0.115(10%) 


0.13 


MWDA [28] 


0.548(< 1%) 


0.474(4%) 


0.156(49%) 


0.10 


RY DFT [29, 30] 


0.60(10%) 


0.511(3%) 


0.174(69%) 


0.06 


Simulation [25,26] 


0.545 


0.493 


0.104 


~ 0.13 




FIG. 1: Density p g (r) of an amorphous structure calculated from Eq.(2.5) using the data of g(R) 
found from molecular dynamic simulation of granular particles subjected to a sequence of vertical 
tapes for a = 150 (strong localization condition) and a = 15 (weak localization condition) 



19 



8 - 

g(r) for isotropic liquid at T)= 0.576 




FIG. 2: Comparison of pair correlation function g(R) of an amorphous structure and homogeneous 
liquid at the same packing fraction 77 = 0.576 
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FIG. 3: Direct pair correlation function c^°\r) as a function of distance r at packing fraction 
r] = 0.50 and 0.55 found from the integral equation theory. The inset shows at the magnified scale 
the value for r > 1 
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FIG. 4: Density derivatives of direct pair correlation function c^°\r) at rj = 0.50 and 0.55 found 
from the integral equation theory. The inset shows at magnified scale the value for r > 1. 
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FIG. 6: Values of harmonic coefficients c^ G \ m (k) for first two sets of reciprocal lattice vectors of a 
f.c.c. crystal for r\ = 0.55 and a = 170 
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FIG. 7: Comparison of order parameters of amorphous structure for a = 50 and 15 and of a f.c.c. 
crystal for a = 50 at r] = 0.596 
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FIG. 8: Values of c^{r) function of distance r at q = 7.53 at which p q is maximum. Inset 
magnifies the value of (r) for r > 1 
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FIG. 9: Free energy difference /S.A = A[po\ — A\[pq\ as a function of localization parameter a at 
rf = 0.576 and 0.596 for amorphous structures. The inset shows at magnified scale the energy 
barrier that separates the minimum of amorphous structure from that of homogeneous fluid. 
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